rm(list = ls())
library(tidyverse)
library(broom)
library(haven)
library(polycor)
library(wCorr)
library(quanteda)
library(quanteda.textplots)
require(glmnet)
require(quanteda.textstats)
library(ggcorrplot)
library(tidymodels)

set.seed(221186)

weightedCorr_na <- function(x,y,weights,...){
  
  x_missing <- which(is.na(x))
  y_missing <- which(is.na(y))
  missing <- unique(c(x_missing, y_missing))
  
  weights <- weights[-missing]
  
  weightedCorr(x[-missing], y[-missing], weights = weights, ...)  
}
mae <- function(x, weights) weighted.mean(abs(x - mean(x, na.rm = T)), na.rm = TRUE, w = weights)

n_boot_reps <- 500

load("../working/survey.Rdata")

## Combine wave 1 respondents with "top up" respondents from wave 2

vars <- c("trans_rights", "offensive_speech", "minimum_wage", "zero_hours", "unemployment_support", "high_tax")

get_effects_all <- function(split, wave = 2, matched_pairs = FALSE){
  cat(".")
  vars <- c("trans_rights", "offensive_speech", "minimum_wage", "zero_hours", "unemployment_support", "high_tax")
  
  if(wave == 2){
    
    for(i in 1:length(vars)) split[[vars[i]]] <- split[[paste0(vars,"_w2")[i]]]

  }
  
  split %>%
  group_by(treat) %>%
  summarise(across(all_of(vars), list(trans_rights = ~weightedCorr_na(.,get(vars[1]), 
                                                              method = "Polychoric", weights = weight),
                              offensive_speech = ~weightedCorr_na(.,get(vars[2]), 
                                                              method = "Polychoric", weights = weight),
                              minimum_wage = ~weightedCorr_na(.,get(vars[3]), 
                                                              method = "Polychoric", weights = weight),
                              zero_hours = ~weightedCorr_na(.,get(vars[4]), 
                                                              method = "Polychoric", weights = weight),
                              unemployment_support = ~weightedCorr_na(.,get(vars[5]), 
                                                              method = "Polychoric", weights = weight),
                              high_tax = ~weightedCorr_na(.,get(vars[6]), 
                                                              method = "Polychoric", weights = weight)), 
                   .names = "{.col}-{.fn}")) %>%
  pivot_longer(-1) %>%
  pivot_wider(names_from = treat, values_from = value) %>%
  mutate(Diff = Treatment - Control) %>%
  filter(Control != 1) %>%
  filter(name %in% apply(t(combn(vars,2)),1,function(x) paste0(x, collapse = "-"))) %>% # Kludge to remove the duplicate correlation pairs
  #separate(name, c("var1", "var2"), "-") %>%
  mutate(name = tools::toTitleCase(gsub("_"," ", name)),
         name = tools::toTitleCase(gsub("-","/", name)),) %>%
    {if(matched_pairs) filter(.,name %in% c("Trans Rights/Offensive Speech", "Minimum Wage/Zero Hours", "Unemployment Support/High Tax")) else .} %>%
  mutate(Diff_mean = sum(Diff)/length(Diff))
  
  
  
}

get_intervals <- function(split, wave = 1, matched_pairs = FALSE){
  
  ests <- get_effects_all(split, wave, matched_pairs = matched_pairs)
  boot_out <- tibble(fx = lapply(1:n_boot_reps, function(x) get_effects_all(split[sample(1:nrow(split), replace = T),], wave = wave, matched_pairs = matched_pairs)))
  
  ints <- boot_out %>% unnest(fx) %>%
    group_by(name) %>%
    summarise(high_treat = quantile(Treatment, 0.975),
              high_control = quantile(Control, 0.975),
              lo_treat = quantile(Treatment, 0.025),
              lo_control = quantile(Control, 0.025),
              high = quantile(Diff, 0.975),
              low = quantile(Diff, 0.025))
  
  average_ests <- boot_out %>% mutate(id = 1:length(fx)) %>% unnest(fx) %>% group_by(id) %>% 
    summarise(Diff_mean = unique(Diff_mean)) %>% ungroup() %>% 
    summarise(Diff_low = quantile(Diff_mean, 0.025),
              Diff_high = quantile(Diff_mean, 0.975)
              )
  
  ests$Diff_low <- average_ests$Diff_low
  ests$Diff_high <- average_ests$Diff_high
  
  ests <- full_join(ests, ints, by = "name")
  
  ests$est <- ests$Diff
  
  ests %>% select(-c(Diff))

}


## Estimate main effects

constraint_effects_w1 <- get_intervals(just_w1, 1, matched_pairs = FALSE)
constraint_effects_panel <- get_intervals(just_panel, 2, matched_pairs = FALSE)
constraint_effects_topup <- get_intervals(just_topup, 2, matched_pairs = FALSE)
constraint_effects_combined <- get_intervals(just_w1_topup_combined,1, matched_pairs = FALSE)

## Estimate paired-effects

constraint_effects_w1_pair <- get_intervals(just_w1, 1, matched_pairs = TRUE)
constraint_effects_panel_pair <- get_intervals(just_panel, 2, matched_pairs = TRUE)
constraint_effects_topup_pair <- get_intervals(just_topup, 2, matched_pairs = TRUE)
constraint_effects_combined_pair <- get_intervals(just_w1_topup_combined,1, matched_pairs = TRUE)

## Estimate attrition-weighted main effects

constraint_effects_w1_attrition <- just_w1 %>% 
  mutate(weight = weight * weight_attrition_in_wave) %>% 
  get_intervals(1, matched_pairs = FALSE)

constraint_effects_panel_attrition <- just_panel %>% 
  mutate(weight = weight * weight_attrition_panel) %>% 
  get_intervals(2, matched_pairs = FALSE)

constraint_effects_topup_attrition <- just_topup %>% 
  mutate(weight = weight * weight_attrition_in_wave) %>% 
  get_intervals(2, matched_pairs = FALSE)

## Plot main effects

constraint_effects_w1$wave <- "Sample One, Wave One"
constraint_effects_panel$wave <- "Sample One, Wave Two"
constraint_effects_topup$wave <- "Sample Two"
constraint_effects_combined$wave <- "Combined Sample"

average_treat_effects <- data.frame(wave = c("Sample One, Wave One", "Sample Two", "Sample One, Wave Two"),
                                    est = c(unique(constraint_effects_w1$Diff_mean), 
                                            unique(constraint_effects_topup$Diff_mean), 
                                            unique(constraint_effects_panel$Diff_mean)),
                                    low = c(unique(constraint_effects_w1$Diff_low), 
                                            unique(constraint_effects_topup$Diff_low), 
                                            unique(constraint_effects_panel$Diff_low)),
                                    high = c(unique(constraint_effects_w1$Diff_high), 
                                            unique(constraint_effects_topup$Diff_high),
                                            unique(constraint_effects_panel$Diff_high)),
                                    name = constraint_effects_topup$name[1]) %>% 
  mutate(wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two", "Sample Two")))

bind_rows(constraint_effects_w1, constraint_effects_panel, constraint_effects_topup) %>%
  mutate(name = factor(name, levels = constraint_effects_w1$name[order(constraint_effects_w1$est)]),
         wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two", "Sample Two"))) %>%
  ggplot(aes(x = est, xmin = low, xmax = high, y = name, col = wave, fill = wave)) + 
  geom_pointrange(position = position_dodge(width = .2)) + 
  facet_wrap(~wave, ncol = 4) + 
  geom_rect(data = average_treat_effects, ymax = 100, ymin = -100, alpha = .4, lwd = 0.0001) + 
  geom_vline(aes(xintercept = est, col = wave), data = average_treat_effects) + 
  theme_bw() + 
  xlab("Effect of reason-giving on constraint") + 
  ylab("") + 
  theme(#panel.background = element_rect(fill = "transparent"),
        #plot.background = element_rect(fill = "transparent", color = NA),
        #legend.background = element_rect(fill = "transparent", color = NA),
        #legend.box.background = element_rect(fill = "transparent", color = NA),
        legend.key = element_blank(),
        legend.position = "bottom") +
  geom_vline(xintercept = 0, linetype = 2) + 
  geom_hline(yintercept = c("Minimum Wage/Zero Hours", "Trans Rights/Offensive Speech", "Unemployment Support/High Tax"), alpha = .2, size = 5) + 
  scale_color_manual("", values = c("black", "black", "black", "black")) + 
  scale_fill_manual("", values = c("black", "black", "black", "black")) + 
  guides(fill = "none", color = "none") 

ggsave("../out/outcomes/constraint.pdf", width = 10, height = 5)
ggsave("../out/outcomes/constraint.png", dpi = 600, width = 10, height = 5)

bind_rows(constraint_effects_w1, constraint_effects_panel, constraint_effects_topup) %>%
  select(name, Control, Treatment, high_treat, high_control, lo_treat, lo_control, wave) %>%
  rename(est_control = Control,
         est_treat = Treatment) %>%
  pivot_longer(cols = -c(name, wave),
               names_to = c(".value", "Var"),
               names_sep = "_") %>%
  mutate(name = factor(name, levels = constraint_effects_w1$name[order(constraint_effects_w1$est)]),
         wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two", "Sample Two")),
         Var = case_when(Var == "control" ~ "Control",
                         Var == "treat" ~ "Treatment")) %>%
  ggplot(aes(x = est, xmin = lo, xmax = high, y = name, col = Var, fill = Var)) + 
  geom_pointrange(position = position_dodge(width = .2)) + 
  facet_wrap(~wave, ncol = 4) + 
  theme_bw() + 
  xlab("Correlation in issue positions") + 
  ylab("") + 
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill = "transparent", color = NA),
        legend.box.background = element_rect(fill = "transparent", color = NA),
        legend.key = element_blank(),
        legend.position = "bottom") +
  geom_vline(xintercept = 0, linetype = 2)  +
  geom_hline(yintercept = c("Minimum Wage/Zero Hours", "Trans Rights/Offensive Speech", "Unemployment Support/High Tax"), alpha = .2, size = 5) + 
  scale_color_manual("", values = c("black", "gray")) + 
  scale_fill_manual("", values = c("black", "gray")) 

ggsave("../out/outcomes/constraint_levels.pdf", width = 15, height = 5)

## Plot attrition-weighted main effects

constraint_effects_w1_attrition$wave <- "Sample One, Wave One"
constraint_effects_panel_attrition$wave <- "Sample One, Wave Two"
constraint_effects_topup_attrition$wave <- "Sample Two"

average_treat_effects <- data.frame(wave = c("Sample One, Wave One", "Sample Two", "Sample One, Wave Two"),
                                    est = c(unique(constraint_effects_w1_attrition$Diff_mean), 
                                            unique(constraint_effects_topup_attrition$Diff_mean), 
                                            unique(constraint_effects_panel_attrition$Diff_mean)),
                                    low = c(unique(constraint_effects_w1_attrition$Diff_low), 
                                            unique(constraint_effects_topup_attrition$Diff_low), 
                                            unique(constraint_effects_panel_attrition$Diff_low)),
                                    high = c(unique(constraint_effects_w1_attrition$Diff_high), 
                                             unique(constraint_effects_topup_attrition$Diff_high),
                                             unique(constraint_effects_panel_attrition$Diff_high)),
                                    name = constraint_effects_topup_attrition$name[1]) %>% 
  mutate(wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two", "Sample Two")))

bind_rows(constraint_effects_w1_attrition, constraint_effects_panel_attrition, constraint_effects_topup_attrition) %>%
  mutate(name = factor(name, levels = constraint_effects_w1_attrition$name[order(constraint_effects_w1_attrition$est)]),
         wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two", "Sample Two"))) %>%
  ggplot(aes(x = est, xmin = low, xmax = high, y = name, col = wave, fill = wave)) + 
  geom_pointrange(position = position_dodge(width = .2)) + 
  facet_wrap(~wave, ncol = 4) + 
  geom_rect(data = average_treat_effects, ymax = 100, ymin = -100, alpha = .4, lwd = 0.0001) + 
  geom_vline(aes(xintercept = est, col = wave), data = average_treat_effects) + 
  theme_bw() + 
  xlab("Effect of reason-giving on constraint") + 
  ylab("") + 
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill = "transparent", color = NA),
        legend.box.background = element_rect(fill = "transparent", color = NA),
        legend.key = element_blank(),
        legend.position = "bottom") +
  geom_vline(xintercept = 0, linetype = 2) + 
  geom_hline(yintercept = c("Minimum Wage/Zero Hours", "Trans Rights/Offensive Speech", "Unemployment Support/High Tax"), alpha = .2, size = 5) + 
  scale_color_manual("", values = c("black", "black", "black", "black")) + 
  scale_fill_manual("", values = c("black", "black", "black", "black")) +  
  guides(fill = "none", color = "none") 

ggsave("../out/outcomes/constraint_attrition.pdf", width = 10, height = 5)

## Plot pair effects

constraint_effects_w1_pair$wave <- "Sample One, Wave One"
constraint_effects_panel_pair$wave <- "Sample One, Wave Two"
constraint_effects_topup_pair$wave <- "Sample Two"
constraint_effects_combined_pair$wave <- "Combined Sample"

average_treat_effects_pair <- data.frame(wave = c("Sample One, Wave One", "Sample Two", "Sample One, Wave Two"),
                                    est = c(unique(constraint_effects_w1_pair$Diff_mean), 
                                            unique(constraint_effects_topup_pair$Diff_mean), 
                                            #unique(constraint_effects_combined$Diff_mean),
                                            unique(constraint_effects_panel_pair$Diff_mean)),
                                    low = c(unique(constraint_effects_w1_pair$Diff_low), 
                                            unique(constraint_effects_topup_pair$Diff_low), 
                                            #unique(constraint_effects_combined$Diff_low),
                                            unique(constraint_effects_panel_pair$Diff_low)),
                                    high = c(unique(constraint_effects_w1_pair$Diff_high), 
                                             unique(constraint_effects_topup_pair$Diff_high),
                                             #unique(constraint_effects_combined$Diff_high),
                                             unique(constraint_effects_panel_pair$Diff_high)),
                                    name = constraint_effects_topup_pair$name[1]) %>% 
  mutate(wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two", "Sample Two")))

bind_rows(constraint_effects_w1_pair, constraint_effects_panel_pair, constraint_effects_topup_pair) %>%
  mutate(name = factor(name, levels = constraint_effects_w1_pair$name[order(constraint_effects_w1_pair$est)]),
         wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two", "Sample Two"))) %>%
  ggplot(aes(x = est, xmin = low, xmax = high, y = name, col = wave, fill = wave)) + 
  geom_pointrange(position = position_dodge(width = .2)) + 
  facet_wrap(~wave, ncol = 4) + 
  geom_rect(data = average_treat_effects_pair, ymax = 100, ymin = -100, alpha = .4, lwd = 0.0001) + 
  geom_vline(aes(xintercept = est, col = wave), data = average_treat_effects_pair) + 
  theme_bw() + 
  xlab("Effect of reason-giving on constraint") + 
  ylab("") + 
  theme(panel.background = element_rect(fill = "transparent"),
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill = "transparent", color = NA),
        legend.box.background = element_rect(fill = "transparent", color = NA),
        legend.key = element_blank(),
        legend.position = "bottom") +
  geom_vline(xintercept = 0, linetype = 2) + 
  scale_color_manual("", values = c("black", "black", "black", "black")) + 
  scale_fill_manual("", values = c("black", "black", "black", "black")) + 
  guides(fill = "none", color = "none") 

ggsave("../out/outcomes/constraint_pairs.pdf", width = 15, height = 5)

constraint_effects_est <- list(constraint_effects_w1 = constraint_effects_w1,
     constraint_effects_panel = constraint_effects_panel,
     constraint_effects_topup = constraint_effects_topup,
     constraint_effects_combined = constraint_effects_combined,
     average_treat_effects = average_treat_effects,
     
     constraint_effects_w1_pair = constraint_effects_w1_pair,
     constraint_effects_panel_pair = constraint_effects_panel_pair,
     constraint_effects_topup_pair = constraint_effects_topup_pair,
     constraint_effects_combined_pair = constraint_effects_combined_pair,
     average_treat_effects_pair = average_treat_effects_pair)

save(constraint_effects_est, file = "../working/constraint_effects_est.Rdata")

## Heterogeneous effects

# Combine wave 1 responses with top-up responses in wave 2

just_topup <- just_topup %>% 
  mutate(high_tax = high_tax_w2,
         zero_hours = zero_hours_w2,
         unemployment_support = unemployment_support_w2,
         trans_rights = trans_rights_w2,
         offensive_speech = offensive_speech_w2,
         minimum_wage = minimum_wage_w2)

tmp <- bind_rows(just_topup, just_w1)

#constraint_effects_attention <- lapply(c("High", "Medium", "Low"), function(x) get_intervals(tmp[tmp$attention_cat == x,], 1))
#names(constraint_effects_attention) <- c("High", "Medium", "Low")

constraint_effects_attention <- lapply(c("High", "Low"), function(x) get_intervals(tmp[tmp$attention_cat2 == x,], 1))
names(constraint_effects_attention) <- c("High", "Low")

constraint_effects_gender <- lapply(c("Male", "Female"), function(x) get_intervals(tmp[tmp$gender == x,], 1))
names(constraint_effects_gender) <- c("Male", "Female")

constraint_effects_age <- lapply(sort(unique(tmp$age_cat)), function(x) get_intervals(tmp[tmp$age_cat == x,], 1))
names(constraint_effects_age) <- sort(unique(tmp$age_cat))

constraint_effects_educ <- lapply(sort(unique(tmp$education.x)), function(x) get_intervals(tmp[tmp$education.x == x,], 1))
names(constraint_effects_educ) <- sort(unique(tmp$education.x))

constraint_effects_eu_vote <- lapply(sort(unique(tmp$eu_vote)), function(x) get_intervals(tmp[tmp$eu_vote == x,], 1))
names(constraint_effects_eu_vote) <- sort(unique(tmp$eu_vote))

constraint_effects_vote19 <- lapply(sort(unique(tmp$vote19_simple)), function(x) get_intervals(tmp[tmp$vote19_simple == x,], 1))
names(constraint_effects_vote19) <- sort(unique(tmp$vote19_simple))

constraint_het_effects <- list(attention = constraint_effects_attention,
                                 gender = constraint_effects_gender,
                                 age = constraint_effects_age,
                                 education = constraint_effects_educ,
                                 eu_vote = constraint_effects_eu_vote,
                                 vote19 = constraint_effects_vote19)

save(constraint_het_effects, file = "../working/het_effects_constraint.Rdata")

# Heterogeneous effects by treatment duration

get_effects_all_min_thresh <- function(split, wave = 2, matched_pairs = FALSE, threshold = 30){
  cat(".")
  vars <- c("trans_rights", "offensive_speech", "minimum_wage", "zero_hours", "unemployment_support", "high_tax")
  
  if(wave == 2){
    
    for(i in 1:length(vars)) split[[vars[i]]] <- split[[paste0(vars,"_w2")[i]]]
    
  }
  
  split %>%
    mutate(trans_rights = ifelse(duration_trans_text < threshold & treat == "Treatment", NA, trans_rights),
           offensive_speech = ifelse(duration_speech_text < threshold & treat == "Treatment", NA, offensive_speech),
           minimum_wage = ifelse(duration_mw_text < threshold & treat == "Treatment", NA, minimum_wage),
           zero_hours = ifelse(duration_zero_text < threshold & treat == "Treatment", NA, zero_hours),
           unemployment_support = ifelse(duration_unemp_text < threshold & treat == "Treatment", NA, unemployment_support),
           high_tax = ifelse(duration_tax_text < threshold & treat == "Treatment", NA, high_tax)) %>%
    group_by(treat) %>%
    summarise(across(all_of(vars), list(trans_rights = ~weightedCorr_na(.,get(vars[1]), 
                                                                        method = "Polychoric", weights = weight),
                                        offensive_speech = ~weightedCorr_na(.,get(vars[2]), 
                                                                            method = "Polychoric", weights = weight),
                                        minimum_wage = ~weightedCorr_na(.,get(vars[3]), 
                                                                        method = "Polychoric", weights = weight),
                                        zero_hours = ~weightedCorr_na(.,get(vars[4]), 
                                                                      method = "Polychoric", weights = weight),
                                        unemployment_support = ~weightedCorr_na(.,get(vars[5]), 
                                                                                method = "Polychoric", weights = weight),
                                        high_tax = ~weightedCorr_na(.,get(vars[6]), 
                                                                    method = "Polychoric", weights = weight)), 
                     .names = "{.col}-{.fn}")) %>%
    pivot_longer(-1) %>%
    pivot_wider(names_from = treat, values_from = value) %>%
    mutate(Diff = Treatment - Control) %>%
    filter(Control != 1) %>%
    filter(name %in% apply(t(combn(vars,2)),1,function(x) paste0(x, collapse = "-"))) %>% # Kludge to remove the duplicate correlation pairs
    #separate(name, c("var1", "var2"), "-") %>%
    mutate(name = tools::toTitleCase(gsub("_"," ", name)),
           name = tools::toTitleCase(gsub("-","/", name)),) %>%
    {if(matched_pairs) filter(.,name %in% c("Trans Rights/Offensive Speech", "Minimum Wage/Zero Hours", "Unemployment Support/High Tax")) else .} %>%
    mutate(Diff_mean = sum(Diff)/length(Diff))
  
  
  
}

get_intervals_thresh <- function(split, wave = 1, matched_pairs = FALSE, threshold = 30){
  
  ests <- get_effects_all_min_thresh(split, wave, matched_pairs = matched_pairs, threshold = threshold)
  boot_out <- tibble(fx = lapply(1:n_boot_reps, function(x) get_effects_all_min_thresh(split[sample(1:nrow(split), replace = T),], wave = wave, matched_pairs = matched_pairs)))
  
  ints <- boot_out %>% unnest(fx) %>%
    group_by(name) %>%
    summarise(high_treat = quantile(Treatment, 0.975),
              high_control = quantile(Control, 0.975),
              lo_treat = quantile(Treatment, 0.025),
              lo_control = quantile(Control, 0.025),
              high = quantile(Diff, 0.975),
              low = quantile(Diff, 0.025))
  
  average_ests <- boot_out %>% mutate(id = 1:length(fx)) %>% unnest(fx) %>% group_by(id) %>% 
    summarise(Diff_mean = unique(Diff_mean)) %>% ungroup() %>% 
    summarise(Diff_low = quantile(Diff_mean, 0.025),
              Diff_high = quantile(Diff_mean, 0.975)
    )
  
  ests$Diff_low <- average_ests$Diff_low
  ests$Diff_high <- average_ests$Diff_high
  
  ests <- full_join(ests, ints, by = "name")
  
  ests$est <- ests$Diff
  
  ests %>% select(-c(Diff))
  
}


constraint_effects_w1_thresh <- get_intervals_thresh(just_w1, 1, matched_pairs = FALSE, threshold = 30)
constraint_effects_panel_thresh <- get_intervals_thresh(just_panel, 2, matched_pairs = FALSE, threshold = 30)
constraint_effects_combined_thresh <- get_intervals_thresh(just_w1_topup_combined,1, matched_pairs = FALSE, threshold = 30)

## Plot pair effects

constraint_effects_w1_thresh$wave <- "Sample One, Wave One"
constraint_effects_panel_thresh$wave <- "Sample One, Wave Two"
constraint_effects_combined_thresh$wave <- "Combined Sample"

constraint_average_treat_effects_thresh <- data.frame(wave = c("Sample One, Wave One", "Sample One, Wave Two"),
                                         est = c(unique(constraint_effects_w1_thresh$Diff_mean), 
                                                 unique(constraint_effects_panel_thresh$Diff_mean)),
                                         low = c(unique(constraint_effects_w1_thresh$Diff_low), 
                                                 unique(constraint_effects_panel_thresh$Diff_low)),
                                         high = c(unique(constraint_effects_w1_thresh$Diff_high), 
                                                  unique(constraint_effects_panel_thresh$Diff_high))) %>% 
  mutate(wave = factor(wave, levels = c("Sample One, Wave One", "Sample One, Wave Two")))

save(constraint_effects_w1_thresh, constraint_effects_panel_thresh, constraint_average_treat_effects_thresh, file = "../working/dur_effects_constraint.Rdata")
